const N: usize = 10000;
static mut PRIMES: Vec<usize> = Vec::new();
static mut NOT_PRIME: [bool; N + 1] = [false; N + 1];
static mut INIT: bool = false;

pub fn init_sieve() {
    unsafe {
        NOT_PRIME[0] = true;
        NOT_PRIME[1] = true;
        for i in 2..=N {
            if !NOT_PRIME[i] {
                PRIMES.push(i);
            }
            for pri_j in PRIMES.iter() {
                if i * pri_j > N {
                    break;
                }
                NOT_PRIME[i * pri_j] = true;
                if i % pri_j == 0 {
                    break;
                }
            }
        }
    }
}

const fn modmul(a: u64, b: u64, modulo: u64, modulo_inv: u64) -> u64 {
    let a = a as u128 * b as u128;
    let (t, f) = ((a >> 64) as u64).overflowing_sub(
        (((a as u64).wrapping_mul(modulo_inv) as u128).wrapping_mul(modulo as u128) >> 64) as u64,
    );
    t.wrapping_add(modulo * f as u64)
}

const fn modpow(x: u64, mut exp: u64, m: u64, modulo_inv: u64, r: u64) -> u64 {
    let (mut res, mut val) = (r, x);
    while exp > 0 {
        if exp & 1 != 0 {
            res = modmul(res, val, m, modulo_inv);
        }
        val = modmul(val, val, m, modulo_inv);
        exp >>= 1;
    }
    res
}

fn miller_rabin(p: u64) -> bool {
    if p == 1 || p & 1 == 0 {
        return p == 2;
    }

    let r = ((1u128 << 64) % p as u128) as u64;
    let r2 = ((p as u128).wrapping_neg() % p as u128) as u64;
    let modulo_inv = {
        let mut inv = p;
        while p.wrapping_mul(inv) != 1 {
            inv = inv.wrapping_mul(2u64.wrapping_sub(p.wrapping_mul(inv)));
        }
        inv
    };

    let s = (p - 1).trailing_zeros();
    let t = (p - 1) >> s;

    [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
        .iter()
        .all(|&a| {
            let a = modmul(a, r2, p, modulo_inv);
            let at = modpow(a, t, p, modulo_inv, r);
            if at == r || at == p - r {
                return true;
            }

            (1..s)
                .scan(at, |at, _| {
                    *at = modmul(*at, *at, p, modulo_inv);
                    Some(*at)
                })
                .any(|at| at == p - r)
        })
}

pub fn is_prime(n: u64) -> bool {
    unsafe {
        if n < N as u64 {
            return !NOT_PRIME[n as usize];
        }
    }
    miller_rabin(n)
}

// 从1开始逆时针螺旋着摆放自然数，我们可以构造出一个边长为7的螺旋数阵。
// 可以发现在所有对角线上一共有 8 个素数，比例达到 8/13≈628/13≈62。
// 如果不断重复这个过程，求对角线上素数的比例第一次低于某个比例时，
// 螺旋数阵的边长是多少，此时有多少个素数？
pub fn min_edge_prime_num(number: u32) -> String {
    unsafe {
        if !INIT {
            init_sieve();
            INIT = true;
        }
    }

    let rate = number as u64;
    let mut i = 1;
    let mut p = 3;
    while (4 * i + 1) * rate <= p * 100 {
        i += 1;
        let rt = (2 * i + 1) * (2 * i + 1) - 6 * i;
        p += is_prime(rt) as u64 + is_prime(rt + (i << 1)) as u64 + is_prime(rt + (i << 2)) as u64;
    }
    format!("{},{}", 2 * i + 1, p)
}
